function [smooth_st,mat_smooth_counterfactual,countmat,mat_countmat_counterfactual,euokmat]=lmj_counterfactual_multiple(funcmod,par1,par2mat,Y,trainvec,into)
% (1) Call LMJ_KFITERCENTRAL using par1

[~,etamat,smooth_st,~,~,~,countmat,~,countszero]=lmj_kfiltercentral(par1,funcmod,Y,trainvec,into.solveopt,into.addsol); 
[T, ns, nx] = size(countmat); %number of shocks in total; 

% check smooth_st and countmat is correctly generated
check1 = nan(nx, 1); 
for ii = 1:nx; 
    check1(ii) = max(abs(sum(countmat(:, ii, :), 3) - smooth_st(:, ii))); 
end
clear ii; 

if ~all(check1<1e-8); 
    error('variable countmat and/or smooth_st not defined corrected')
end

% Check feeding ETAMAT recovers observables and states 
disp('Checking subroutine...'); 
[ycounter,stcounter,countobs_check,countmat_check]=lmj_counterfactual_sub(funcmod,par1,into.solveopt,into.addsol,countszero,etamat);
maxdif=comparemat(Y,ycounter); 
if maxdif > 1e-8; 
    error('COUNTERFACTUAL_SUB is not working correctly on Y') 
end 
maxdif=comparemat(smooth_st,stcounter); 
if maxdif > 1e-8; 
    error('COUNTERFACTUAL_SUB is not working correctly on S') 
end 

nselect=length(into.state_pos); 
numPar = size(par2mat, 2); 
mat_smooth_counterfactual  = nan(T, nselect, numPar); 
mat_countmat_counterfactual= nan(T, nselect, nx,  numPar); 
euokmat = nan(numPar, 1); 
for jj = 1:numPar
    par2 = par2mat(:, jj); 
    disp(' ')
    disp('___________________________________')
    disp(['Starting Parameter ', num2str(jj)])
    if into.flag_mixfreq==0
        
        disp('FIX ME'); 
        error('FIX ME'); 
        % (2) solve model using par2
        [GG2, RR2, CONS2, eu2, SDX2, ZZ2]=feval(funcmod,par2,varargin{:});
        if ~isequal(eu2, [1;1]);
            disp(['Cannot find a solution for par ', num2str(jj)])
            euokmat(jj) = 0;
            continue
        else
            euokmat(jj) = 1;
            disp(['Solution found for par ', num2str(jj)])
        end
        smooth_counter = nan(T, ns);
        smooth_counter(1, :) = sum(countszero, 2); % initialization
        for tt = 2:T;
            smooth_counter(tt, :) = (GG2*smooth_counter(tt-1, :)'+ RR2*etamat(tt, :)')';
        end
        mat_smooth_counter(:, :, jj) = smooth_counter;
        
        counter_countmat = nan(T, ns, nx);
        for ii = 1:nx;
            tempeta = zeros(nx, T);
            tempeta(ii, :) = etamat(:, ii)';
            stsim = nan(ns, T);
            stsim(:, 1) = countszero(:, ii);
            
            
            for tt = 2:T;
                stsim(:, tt) = GG2*stsim(:, tt-1) + RR2*tempeta(:, tt);
            end
            counter_countmat(:, :, ii) = stsim(:, :)';
        end
        mat_counter_countmat(:, :, :, jj) = counter_countmat;
        % check smooth_counter and counter_countmat is correctly generated
        check2 = nan(nx, 1);
        for ii = 1:nx;
            check2(ii) = max(abs(sum(counter_countmat(:, ii, :), 3) - smooth_counter(:, ii)));
        end
        if ~all(check2<1e-8);
            error('variable countmat and/or smooth_st not defined corrected')
        end
    else 
        [~,stcounter,~,countmatS]=lmj_counterfactual_sub(funcmod,par2,into.solveopt,into.addsol,countszero,etamat);
        if ~isempty(ycounter);
            mat_smooth_counterfactual(:,:,jj)=stcounter(:,into.state_pos);
            mat_countmat_counterfactual(:,:,:,jj)=countmatS(:,into.state_pos,:);
            euokmat(jj)=1; 
        end        
    end
end